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The exact computation of orbits of discrete dynamical systems on the interval is considered. There- 
fore, a multiple-precision floating point approach based on error analysis is chosen and a general 
algorithm is presented. The correctness of the algorithm is shown and the computational complexity 
is analyzed. As a main result, the computational complexity measure considered here is related to 
the Ljapunow exponent of the dynamical system under consideration. 

1 Introduction 

Consider a discrete dynamical system {D,f) on some compact interval D C M, called the phase space, 
given by a function f : D ^ D, a. recursion relation x^+i = f{xn) and an initial value xq € D. The sequence 
{xn)n of iterates is called the orbit of the dynamical system in phase space con^esponding to the initial 
value xq. If such a dynamical system is implemented, that is a computer program is written for calculating 
a finite initial segment of the orbit for given aq, care has to be taken in choosing the appropriate data 
structure for representing real numbers. Traditionally, IEEE 754 double floating point numbers ifTOl 
are used. However, if the dynamical system shows chaotic behavior, a problem arises. The finite and 
constant length of the mantissa of a double variable causes round off errors which are magnified after 
each iteration step. Only after a few iterations, the error is so big that the computed values are actually 
useless fV^. To put things right, a rigorous method for computations with real numbers has to be used. 
In 1 2 1, this issue is addressed for the logistic map which is also consiedred as a starting point in the next 
section. There, the exact real arithmetic in the form of centered intervals with bounded error terms is 
used as described in [31. However, the notation used in IHO is an algebraic one based on arbitrary large 
integers. On the other hand, the aim of the present paper is to keep the notation as close as possible to 
the standard in scientific computing but being precise in the sense of exact real arithmetic. This has as 
a consequence first that the basic data type is not an integer as considered in [2, 3], but a floating point 
number with a definite mantissa length. Second, the type of error considered here is the relative error as 
is standard in floating point arithmetic - in contrast to the absolute error considered in In practice, 

a multiple precision floating point library providing floating point numbers with arbitrary high mantissa 
length have to be used. In the following, it is analyzed how the needed mantissa length behaves in 
multiple-precision computations of iterates of discrete dynamical systems. The mantissa length needed 
for floating point numbers such that any computed point of the orbit has a specified and guaranteed 
accuracy is examined. Therefore, a precise mathematical framework for floating point computations 
has to be established. The main result shows that the ratio of mantissa length to iteration length in 
the limit of iteration length to infinity is related to the Ljapunow exponent. Comparing, in [2] only 
the logistic map is considered explicitly and the connection to the Ljapunow exponent is not stated. 
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but observed numerically. In the present paper, this connection is shown mathematically for a general 
discrete dynamical system {D,f). This result also gives some advice for economically designing exact 
algorithms simulating one-dimensional discrete dynamical systems. 

2 Roundoff Error, Error Propagation and Dynamic Behavior 

In this section, the discrete dynamical system {D,f^) with D = [0, 1] and f^.D^D, f^{x) := /ix(l —x) 
for some control parameter n € (0,4] is investigated. In the literature, the recursion relation = 
fnixn) is called the logistic equation [5|. When implementing the logistic equation on a real computer 
and demanding to obtain exact values for the orbit (x„)„, the analysis of roundoff errors and of error 
propagation requires some care. This is due to the fact that for some values of jU the dynamics is highly 
chaotic and therefore inaccuracies are magnified exponentially in time ||6l|3- 

In the following, for a given initial value xq, the true orbit is denoted by (a:„)„, whereas the really 
computed orbit, suffering from roundoff errors and error propagation, is denoted by (xn),,- Note that 
even xq may differ form xq since the conversion to a floating point number may cause the very first 
roundoff error. One goal of this section is to give a rigorous estimation of the total error in dependence 
of the iteration step n. 

Calculating the orbit (Jc„)„, two types of error are present. First, error propagation due to the iteration 
scheme and second the roundoff error caused by the calculation of f^. Now, let x^ for some n € N be 
given. Then the true error after one iteration step is Xn+i —x„+\. Since in reality not f^{x„) is calculated 
but some erroneous approximation (xn), the true error can be written as x„^\ —Xn+i = {xn) — {x„). 
Hence, the true error can be written as a sum 



of two terms. The first term describes solely the error propagation while the second term gives exactly 
the newly produced error due to the approximate calculation of . 

To handle the exact values of both errors computationally, interval arithmetic can be used HI . Interval 
arithmetic can be seen in the setting here as a special case of the computational model of TTE |fT6l , which 
gives a precise notion for describing computations over the real numbers. Another strongly related model, 
which in some sense reflects the situation here more adequate is the Feasible Real RAM model [4]. For 
the sake of simplicity however, an interval setting is used here. For any time step n, let the phase point 
Xn together with its error be represented by two floating point numbers xl^ and x" {x" > x^) with given 
mantissa length m„ forming an interval [x'f^,x'^]. The interval is an enclosure of the real value x„, that is 
Xn € [xj,,x"]. It is straightforward to transform the interval to a floating point value x„ of mantissa length 
m„ by setting 



where gl{.) performs the rounding to nearest floating point number. The absolute error en '■= |x„ — x„| of 
Xn can be estimated via the interval length dn '■= x" — x', by 



Xn+ l-Xn+l = ifn (X„ ) - //i (X„ ) ) + (/^ (x„ ) - /;j (x„ ) ) 



(1) 




(2) 



en < -^dn + r, 



n 



(3) 



where r„ is an error introduced by the rounding operation gl{.) in Equation |2] An upper bound on r„ will 
be discussed later, for now it suffices to say that in general it is small compared to dn- 
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For doing an error analysis of the logistic equation analytically, some idealizing assumptions are 
made. First, the value of n is assumed to be given with such a high precision that no interval repre- 
sentation is needed. Second, only the error propagation is considered caused by the initial error due to 
rounding xq to some floating point number of mantissa length m. Third, the value of r„ in Equation [3] is 
neglected. The recursion relation then reads in natural interval extension 

with the interval length d„ given by the recursion relation 

dn+ 1 = J — X,,^ J = /I (x„ — — X,^ + X^X^ ) 

= l^dn 

with the obvious solution dn = pi"dQ. Finally the absolute error €„ of x„ according to Equation [2] can be 
bounded from above by 

en < = ^l^"do. (4) 

Note that the above analysis only holds if the natural interval extension for is derived from the formula 
/xx( 1 — x) . If it is derived from the formula (x) = /i (x — x^ ) or (x) = j — (x — ^ )^ , the mathematical 
analysis is more difficult. However, the problems described in the following also appear, but in some 
different form. 

The aim now is to calculate, for given € N, p G Z and mantissa length m, the orbit up to time 
with relative error 10^''. That is, for (x„)o<h<a' should hold 

e„ = |x„ -Xn\< lO^Pxn < IQ-P. (5) 

The ideal assumptions require the somewhat unreal setting that the mantissa length has to be set to some 
finite, but big enough value m for representing xq and a virtually infinite value nioo for doing the iteration. 
Finally, some upper bound on do is needed. The value of do is given as the roundoff error by representing 
Xq as a floating point number of mantissa length m. For that, the well known estimate 

do < 2-'"+^xo < 2-'"+! (6) 

exists. Combining (01) and Q gives as a sufficient condition 

pi" <\0-P (7) 

for?2 = 0, ...,A'^. 

The minimal m, fulfilling the precision requirement ^ on the relative error of x„, which depends 
on xo, A'^ and p, is denoted by m,„in{xo,N,p). So, the sufficient condition ^ gives an upper bound on 

{xo,N,p) by 

m^Uxo,N,p) < r7?-ld(10)+A^-max(0,ld(jU))l (8) 

where ld(.) is the logarithm to base 2. At that stage, a central quantity of this work is introduced which 
is a kind of complexity measure. The loss of significance rate o{x,p), which may depend on the initial 
value X = Xq and the precision p is defined by 

I ^ ,. m,„in{x,N,p) 
o[x,p) :=limsup . 

N^oo N 
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This quantity describes the limiting amount of significant mantissa length being lost at each iteration 
step. Significant means here the part of the places being exact. A general treatment of this complexity 
measure is given in the next section. Roughly speaking, \a{xo,p)N + p ■ Id(lO)] is the mantissa length 
for any floating point number needed in an algorithm doing the iteration starting with xq and calculating 
up to xn, if the output should be precise up to p decimal places. Formula [8] gives an upper bound for the 
loss of significance rate by G{x,p) < max(0,ld(/i)). 

It is interesting to see whether the upper bound calculated analytically, which needed strong ideal- 
izations, is in the region of the real value. So, the logistic equation was implemented using a multiple- 
precision interval library. For that purpose, the interval library MPFI lITSl based on the multiple-precision 
floating point number library MPFR |l8l, both written in C, was used. For each control parameter jj. 
ranging from 0.005 to 4 and a step size of 0.005, the orbit for initial condition 0.22 was calculated up 
to N = 2000. For each jj., the minimum mantissa length m,„;„ needed to guarantee e„ < 10^^;c„ for 
n = 0, . . . ,N was searched. Then, Oest '■ = m,nin/N was calculated. The result shows that Oest exceeds the 
analytical bound max(0,ld(;U)) only slightly. So, the above made ideal assumptions seem to be valid. 
In (VT\, the logistic equation was also investigated for p =3.75 using the exact real arithmetic package 
iRRAM based on the Feasible Real RAM model [4]. In the paper, also the precision needed to guar- 
antee the exactness of the first 6 decimal places are reported up to A'^ = 100000. The values are in full 
agreement with the simulation results performed here. 

Hence, for /i > 1, the interval length increases exponentially in time n. This result should be 
interpreted in terms of the dynamical behavior of the logistic equation. So, at this point is worth having 
an analytical look at the behavior of the dynamical system. Despite the fact that these results are well 
known |I9ll21^ they are reviewed here for the sake of self containment. First, the equation possesses in 
the range D = [0, 1] exactly one fixed point = if /X S (0, 1] and exactly two fixed points x" = and 
x'^^^ = 1 - if /I G (1,4]. Since /^(O) = p and fl{x^^^) = 2 - ^u, x° is a stable fixed point (an attractor, 
\f'ii{^)\ < 1) for M ^ (0: 1) and an unstable fixed point (a repeller, |//,(x'')| > 1) for p £ (1,4]. If = 1, 
the only fixed point x" is hyperbolic (l/i'(^°)| = 1) and a bifurcation occurs at that value of the control 
parameter p. If p €^ {\,3), x" becomes unstable and the newly occurring fixed point x^^^ is stable. Finally, 
lim„^oc/^(;c) = jc^^) for > 1 and lim„^oo/^(x) = x" if p < I holds for all x e (0,1). If p e (0,1), this 
is a direct consequence of the contraction mapping principle. If p = I, observe that /i (x) < x holds for 
all X G (0, 1). Hence, any sequence (/"(x))„, x E (0, 1), is strictly decreasing and bounded from below. 
So it converges to the only fixed point x". For the case p G (1,3), the interested reader is referred to the 
literature: [7 1, Proposition 5.3. At jLt = 3 a second bifurcation occurs and for /i > 3 the system goes into a 
region of periodic behavior with period doubling bifurcations. Finally, for some p < 4, chaotic behavior 
is reached. 

This analysis shows that in the parameter range p G (0, 3), the orbit converges to the stable fixed point 
for any initial value xq G (0, 1). Furthermore, there exists some closed interval I QD, which depends 
on p, containing the stable fixed point such that /^(/) C / holds and is a contraction on /. The 
interval computation using a natural interval extension of the recursion formula px{l —x), on the other 
hand, is not very compatible with this picture. While for p G (0, 1), the results are in agreement with 
the dynamical analysis, the calculations for /i G (1,3) are not handled very well by interval arithmetic 
since the interval approach would suggest an exponential divergence of initially nearby orbits which is 
not true in reality. The reason is that the natural interval approach implicitly, due to the dependency 
problem, takes account only of the global behavior of in the form of the global Lipschitz constant 
max{]/j^(x)] : x G D} = p. However, the local Lipschitz constant max{]/j^(x)] : x G [x|,,x,"]} governs the 
real error propagation at time step n and also describes the dynamic behavior. This notion can be made 
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precise and finally leads to a more efficient algorithm for computing orbits. 

Let us return to Equation [T] The true error is the sum of the error propagation (first term) according 
to the iteration and the roundoff error due to the computation of (second term). The first term of 
Equation [T] can be handled using the mean value theorem, \f^{xn) — f^{xn) \ = \ fliiyn) \ • \xn —Xn\ with 
y,! S [xn — e„,f„ + en]. This gives directly the bound 

\f^i{x„) - f^{xn)\ < sup(|/^([i„-e„,x„+f'„])|)e„. 

The second term can be estimated the following way. As discussed in ifTTl . the roundoff error produced 
in calculating can be estimated by 

\Mx)-Ux)\<l.06K2-'"\f^{x)\ 

where K is the number of rounding operations performed in computing and m is the mantissa length 
of X. In the case considered here, K = 4 follows since there are 3 arithmetic operations and the rounding 
of jJ.. It is further crucial to mention that the factor 1.06 is only valid if < 0.1 • 2'" holds so that the 
mantissa length must not be chosen too small. Using the fact that (x) < ^ holds and (x) < x if 
jJ. <l, the unknown value \f^{x) \ can be eliminated. This calculation shows that there exists a recursive 
equation on an upper bound e„ on e„ for all n: 

Cn+i = L{xn,e„)en + \ .06K2-'"E ^{x„) , eo = 2"'" (9) 

withL(x,e) ■.= sup{\f'^{[x — e,x + e])\) and 



Efi{x) 



ix if /I < 1 
f ifAt>l' 

The idea is now not to calculate intervals, but pairs of values x„ and corresponding guaranteed error 
bounds en- The difference to the interval concept is not to compute the errors implicitly, so that only 
global behavior can be taken into account, but to compute them explicitly and independent of the values 
of interest. It should be mentioned that the approach described here is compatible with an interval 
approach using special centered forms, namely mean value forms |[T4l . However, the approach here 
explicitly devises values and errors, describes an automated error analysis, whereas an interval approach 
primarily does not disclose any error. Furthermore, also the iRRAM package permits a more elaborate 
way for computing the iteration, based on a similar algorithm as described above [13|. The rounded 
values Xn are calculated as usual in floating point arithmetic except that multiple-precision floats are used. 
The guaranteed error bounds are also calculated using floating point according to where interval 
arithmetic is used for calculating L. Only standard precision is needed for calculating the error bounds. 
Implementing this improved algorithm using MPFR and MPFI, the setting as given in the interval case 
produces the following result. In the parameter range jJ. G (0,3), the dynamic behavior is reflected very 
well. Furthermore, in the range p. € [3,4], the curve suggests a relation between the loss of significance 
rate and the Ljapunow exponent A(x) for the logistic map (for a curve of the Ljapunow exponent of the 
logistic map see |5|): a(x) = max(0,A(x))/ln(2) for all p G (0,4]. To be complete, the definition of the 
Ljapunow exponent reads 

Definition 2.1. Let {D,f) be a dynamical system, D C M compact and f : D ^ D continuously differen- 
tiable on the interior of D. Then the Ljapunow exponent at x is defined by 

A(x):=limi£ln|/(/(x))| (10) 

if the limit exists. 
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The Ljapunow exponent may depend on x. However, the following properties hold: 

(a) If {D,f) has an invariant measure p, then the limit in Equation [TOl exists p -almost everywhere. 

(b) Furthermore, if p is ergodic then A(x) is p -almost everywhere constant and equal to 



These properties are a direct consequence of the Birkhoff ergodic theorem, see 1,11,1 . Theorem 4.1.2 and 
Corollary 4.1.9. 

3 The General Algorithm and its Complexity 

Let D be a compact real interval and / : D — )• D a self mapping. In the following, / is assumed to be 
continuous on D, continuously differentiable on the interior of D and /' is bounded. Furthermore, / and 
/' are assumed to be computationally feasible. The precise definition of "computationally feasible" is 
given below. 

In this section, a general algorithm for computing the iteration 



is presented. To be more precise, for given N G'N and € Z, this algorithm computes a finite part of 
the orbit, {xn)Q<n<N, exact in the sense that the relative error at each point x„ does not exceed 10^''. 
The correctness of the algorithm and its computational feasibility is shown. Finally, its complexity is 
examined. 

3.1 Syntax, Semantics and the Algorithm 

The set of all computationally accessible real numbers are the floating point numbers of arbitrary man- 
tissa length denoted by M. In the following, by a floating point number any real number is meant which 
can be expressed by normalized scientific notation. Hence, the set M C M of all floating point num- 
bers is countable infinite and therefore a natural basis for standard computability considerations. Let 
jc G M be some floating point number, then x has as an essential property, its mantissa length denoted by 
x.m. Any real number x is represented in an algorithm concerning real computation by a pair [x] € W- 
consisting of a floating point number [x] .fl approximating x and an upper bound on the relative error, 
[xj.err > 0, also being a floating point number. Furthermore, the inequality \[x].fl —x\< [x].err holds. 
The pair [x] is called a finite precision representation of x. Although [x] .err has the property mantissa 
length, it is irrelevant in what follows. So, the mantissa length of [.ij.err can be assumed to be some 
big enough constant value. Analogously, a function / : D ^ D, Z) C M, is called computationally fea- 
sible if a pair [/] exists of a computable (partial) function [/] .// : R — > M approximating f on D and a 
computable (partial) function [f].erf : — )• M giving an upper bound on the absolute error of [/].// in 
the sense \[f].fl{[x].fl) — f{x)\ < [f].erf{[x]). Here, a partial function / : R ^ M is called computable 
if / is computable as a string function over some finite alphabet where the floating point numbers are 
interpreted as finite strings. Finally, computability over integers, computability of functions with mixed 
arguments and computable predicates are defined in a standard way. 
The algorithm with the above described specification reads 




Xn+l=f{Xn), XqGD 



(11) 
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1 Input parameter: XQ,N,p 

2 Initialize mantissa length m^mo 

3 do 

4 Initialize value and error [x] -(^ gl{xo,m) 

5 for« = OtoA^do 

6 If prec{[x],p) is true then 

7 If not printed print n, [x].fl, [x].err 

8 else break 



10 end for 

11 [x].fl.m^[x].fl.m+\ 

12 while prec{[x],p) is false 

To initialize [x], a rounding function : M x N ^ is needed where gl{xo,mo).fl is a floating point 
number of mantissa length mo being the exactly rounded value of xq to some rounding convention. 
Clearly, gl{xo,mo).err is an upper bound on the absolute rounding error, e.g. gl{xQ,mQ).err = ^ulp{xo) 
if the rounding mode is to nearest. The predicate prec : x Z ^ {true, false} is a test whether the 
relative error of [x], \[x].fl —x\/\x\ if ;c 7^ 0, is bounded by 10'^. The semantics reads: If [x] S is a 
finite precision representation of a' € M and prec{[x] , p) = true holds, then | [x] .fl —x\ < 10^'' \x\ follows. 

In the following, some abbreviations are used occasionally. The floating point numbers and functions 
are indicated by a hat: x := [x].fl and / := [/].//. An over-bar indicates an error bound: e := [x\.err and 
erf := [f].erf. Hence, [x] is equivalent to {x,e) and [/] is equivalent to {f,erf). 

Finally a remark on optimization. The algorithm is not optimized in the performance. Otherwise, 
in Line 10 something like m 2m should be used. Here, the aim is to find the minimal m to guarantee 
some given upper bound on the relative error of x„. 

3.2 Feasibility and Correctness 

It is clear, that the rounding function gl is computationally feasible. So lets begin with the predicate prec. 
Proposition 3.1. The computationally feasible formula 



fulfills the above described semantics. 

Proof. Let {x,e) be a finite precision representation of x. So, if e < 10^^|.x:| holds, then also |jc — ;c| < 
IO~^|.x:| holds. If {x — e){x + e) > 0, then |i| — ^ < \x\ holds. Hence, if {x — e){x + e) > and e < 
lO~P{\x\ — e) holds, then also \x — x\ < 10^^|;t:|. Finally, if ^ < j|%-p |^| holds, then also {x — e){x + e) > 
0. 

Formula [12] only uses the accessible floating point values x and e, basic arithmetics and finite tests. 
Hence, this formula is computationally feasible. □ 

Note that the definition of the predicate this way also gives true in the singular case where x = and 
e = and hence x = 0. 

An algorithm for computing / is by assumption possible. To derive an algorithm for computing erf 
on the absolute error, return to Equations [T] and |9l 



9 



M^[/](H) 




(12) 
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Proposition 3.2. Assume that f{x) computes f{x) up to a correctly rounded last bit in mantissa accord- 
ing to rounding convention. Then there exists a constant K > such that the absolute error of f{x) of 
the computation [/] ( [x\ ) is bounded from above by 

L{x,e)-e+ ^_^^_J f{x)\ (13) 

ifK2^'" < 1. Here, L{x,e) := sup(|/'([^ — ^,^ + ^])|) and m is the mantissa length ofx: x.m. 
Furthermore, this bound is computable. 

Proof. Using Equation [U and following the calculations leading to equation |9l |/(x) — f{x)\ < L{x,e) ■ 
e+ \f{x) — f{x)\ follows. According to the assumption on /, |/(x) — f{x)\ < K2^"'\f{x) \ holds, with 
a value K G {1,2} depending on the rounding convention. However, f{x) is unknown, only /(x) is 
accessible. To overcome this, set /(x) — f{x) = 5f{x) with |5| < K2^"\ Since |5| < 1 holds, resolve to 
/(•^) = TTj/(-^)- Hence, 



\fm-fix)\ 



1 + 5 



f^'y—in 

\m\<j-^jm 



follows. Since an upper bound on L{x,e) can be computed using global optimization techniques, e.g. 
with interval arithmetic, the above described bound is computable. □ 

To summarize, the mathematical iteration (ITTI ) is performed in the algorithm by iterating a value x„ 
approximating x„ with an upper bound on its absolute error e„ according to 

Xn+l=f{Xn) XQ=gl{xQ,m) (14) 

e„+i=L[Xn,en)en + Y—j^;p;;;\Xn+i\ = jZ^kt^^^'^^ ^ ' 

where L(x„,^„) is computable upper bound on L{xn,e„) as described in the preceding proposition. This 
is Line 9 in the inner f or-loop of the algorithm which is executed with successively increasing man- 
tissa length m, controlled by the outer do-while-loop. Finally, it has to be shown that this outer loop 
eventually terminates. Therefore, two more propositions are needed. 

Proposition 3.3. Let x ^ be a real number and {[x\m)m>mo ci sequence of finite precision representa- 
tions of X with increasing mantissa length {[x\m)-fl-m > m such that limm^oo{[x]m) -err = holds and 
consequently limm^oo([x]m).// = x. Then limm^oo prec([x]m , p) = true follows for all p G Z. 

Proof. Since X 7^ there exists some M € N such that < ^|x| < |([x]„,).//| and {[x]m) -err < 2(1+10-/') l-^l 
holds for all m > M. Then, prec{[x\m, p) = true holds for all m>M. □ 

The next proposition makes the link to Line 9 in the algorithm. 

Proposition 3.4. LetXn be the n-th element of the orbit ofEquationMl \and {[xn]m)m>mo ^ sequence given 
according to the recursion equations di4D and di5D with increasing mantissa length ([x„],„).fl.m > m. 
Then lmim^c»([x„],„).err = holds and consequently \im,„^oo{[xn]m) -fl = ^n- 

Proof. Let L := sup (/'(D)) and L > L be some computationally accessible value using some global 
optimization technique. Then Equation [T5] leads to e„+i < Le„ + ^^^^„, M where M > sup{[x| : x € D} 

such that |x„| < M holds for all n. Iteration gives e„ < L"eo + j^^^MYJI^^l'^ < j^^^MJJI^qL!^ . 
Hence, for n fixed, limm^oo([x„]„,).err = follows. □ 
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These two propositions finish the correctness proof of the algorithm. They show that, if x„ 7^ for 
n = 0,...,N, the outer loop eventually terminates for any G Z. 

3.3 Computational Complexity 

After having presented the preliminary work, the main issue of the paper is addressed - the computational 
complexity of the presented algorithm. The complexity measure of interest here is the loss of significance 
rate already introduced informally in the last section. Here is the formal definition. 

Definition 3.1. The minimal mantissa length, for which the described algorithm eventually halts is de- 
noted by m,„in{xQ,N ,p), where xq, N and p are the corresponding input parameters. Then, the loss of 
significance rate a:]RnDxZ^>Mis defined by 

^ ,. m,„in{x,N,p) 
a(;ic,p) := limsup — . (16) 

However, to achieve bounds on the loss of significance rate, a technical difficulty has to be circum- 
vented. Therefore, one more assumption on the dynamical system {D,f ), additional to the ones already 
mentioned in the beginning of this section, has to be made. 

Assumption 3.1. The dynamical system {D,f) is assumed to have the properties already mentioned in 
the beginning of this section and additionally ^ D. 

It was already seen in the last subsection that jc,, = makes difficulties such that it cannot be proven 
that the algorithm eventually halts. However, the restriction ^ D is no loss of generality. If all other 
conditions are fulfilled except that D contains zero, consider the following dynamical system {D,f) 
instead. Choose some M > min(D) and set D := {x + M | x G D} as well as f{x) := f{x — M) +M for 
all X £ D. Then {£>,/) fulfills all required properties. Furthermore f'{x) = f'{x — M) holds and therefore 
there is no substantial difference in the complexity analysis of the algorithm between the original system 
and the modified system. 

First, the boundedness of o{x) is shown. 

Proposition 3.5. Let {D,f) be CIS in A^ssuffiptiofi \3. l\ cind fyimin 

{xo,N,p) as in Definition\3l\ Then, for 
given p £lj, there exist soifi€ Ci,C2 ^ 0, dependent of f, such that niinii^lxQ^ N,p) <CiN + C2 holds for 

allN en, xeRnD. 

Proof. According to the requirements made on {D,f), there are some constants L > and M > such 
that e„+i < Len + 1^x2-"' ^ holds for all « € N and all mantissa lengths m. Without loss of generality 
assume L 7^ 1, otherwise set L > 1. Analogous to the treatment in the proof of Proposition 13. 4[ iteration 
gives eN < j^^^MYJ^^^qL" = i^\2-"' ^ ^'^L-\^ ■ Since there exists some B > with B < \x„\ for all 
n, cn/IxnI < eN/B < C2-'"L^+^ follows with C := MK/{B{1 -K2-'^){L- 1)) where mq is the initial 
mantissa length. Line 2 in the algorithm. Then, if C2^'"L'^^^ < ^^^q-,, holds, prec{{xn,e„),p) = true 
for sAXn = 0, . . . ,N . This leads to mmm{xo,N,p) < 

max{0,ld{L))N + max{mo,ld{L) + ld{C) +/? • Id(lO) +ld(l + 10"'')). □ 

Corollary 3.1. Let {D,f) be as in Assumption \3. 1 1 and <j{x,p) the loss of significance rate. Then, for 
given p E Z, there exists some constant C > such that a{x,p) < C holds for all x eM^DD. 

The treatment has now come to a stage that the main statements of this paper can be formulated. 
A lower and an upper bound for the loss of significance rate is given. Furthermore, these bounds are 
strongly related to the Ljapunow exponent A(x) defined in the previous section. 
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Theorem 3.1. Let {D,f) be as in Assumption \3.1\ o{x,p) the loss of significance rate and X{x) the 
Ljapunow exponent of{D,f). Then (7{x,p) > X{x)/ln(2) holds for all x € MPlD, G Z if X{x) exists. 

Proof. First there are two constants B,M >0 such that |x„| > B and \x„\ <M holds for all n. According 
to Equation [T5] and Proposition 13.21 e„^i > |/'(x„)j^„ holds. Iteration gives ejs/ > i^^2-'" ^n=o \f'i^n)\- 
Hence, ^ > j^^^^f^Yln^o \f'i^'')\ follows. A necessary condition for the algorithm to terminate is 
therefore ^2-'"Un=o < TTIo^ which gives m^in{xo,N,p) > Ln^^ ld{\f {x,)\) + p ■ Id(lO) + 

Id(^) + ld(l + 10^''). Following the definitions of the loss of significance rate and the Ljapunow expo- 
nent, a{xo,p) > A(xo)/ln(2) follows. □ 

Before a realistic upper bound on the loss of significance rate can be presented, one more definition 
is needed. 



ln(;t:) ifx>a 
ln(a) ifx<a 



Definition 3.2. Let a > then define a function rja '■ (0,°°) ^ M by 

Furthermore, for any a > define 

Ia{x) :=limsup-£T]„(|/'(/(x))|) 



«— 1 



Proposition 3.6. For all (X> there exists some constant C > such that (x) < C holds for all x € D. 
Furthermore, if the Ljapunow exponent A(x) exists, A(x) < Ao;(x) holds. 

Proof Let L be the Lipschitz constant of / and a > 0. Then for all n£N, ^I^Iq T]a(|/'(/(-^))|) < 
ln(max(a,L)) holds. Hence it follows limsup„^„^^^^Q77a(|/'(/''(x))|) < ln(max(a,L)). The second 
assertion follows from the fact that ln(x) < Tjo;(x) holds for all x > 0, a > 0. □ 



Proposition 3.7. Let x £ D be given. IfX{x) exists, then also the limit 

limXaW=:I(x) (17) 

O!\0 

exists and X{x) > X{x). 

Proof. Since ln(x) < Tja(x) < rip{x) holds for allx > 0, < a < j3, also A(x) < Xa{x) <Xp{x) follows. 
So if a converges in a monotonic decreasing way to 0, a \ 0, the assertion follows. 

□ 

Theorem 3.2. Let {D,f) be as in Assumption \3.1\ a{x,p) the loss of significance rate and X{x) as 
in f lTTl ). Let x G M n D fte given, then for any e > there is some po £ Z such that for all p > po, 
(^(^tP) ^ X{x)/ln{2) +e holds ifX{x) exists. 

The proof is similar to the proof of Theorem 13.11 and can be found in the full version of this article 

mi. 



Christoph Spandl 
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4 Conclusions 

In this paper, two main issues are addressed. First it is shown that a mathematically precise treatment 
of multiple-precision floating point computability is not hard to do. Furthermore this treatment is in a 
manner which is familiar to people working in the field of numerical analysis or scientific computing and 
also for theoretical computer scientists. Furthermore, the formalism does not only allow exact answers 
concerning the existence of a computationally feasible algorithm, but is also allows a treatment of its 
complexity. As a consequence, the described algorithm is formulated not only in an exact and guaranteed 
way, but also enables a motivated reader the real implementation and gives a practical performance 
analysis. 

Second, the results show that the Ljapunow exponent, a central quantity in dynamical systems theory, 
also finds its way into complexity theory, a branch in theoretical computer science. In dynamical systems 
theory, the Ljapunow exponent describes the rate of divergence of initially infinitesimal nearby points. 
For two points having a small but finite initial separation, the Ljapunow exponent has only relevance for 
short time scales [6|. The reason is that due to the boundedness of D, any two different orbits cannot 
separate arbitrarily far away. However, the loss of significance rate shows that the Ljapunow exponent 
has on long time scales not only an asymptotic significance but also a concrete practical one. 
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